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Atmospheric  Structure  Simulation 
An  Autoregressive  Model  For  Smooth 
Geophysical  Power  Spectra  With  Known 

Autocorrelation  Function 


1.  INTRODUCTION 

Atmospheric  fluctuations  in  wind  speed  and  gravity  waves  are  characterized  by  continuous 
power  spectral  density  functions.  For  example,  one  dimensional  wind  speed  PSD's  are  found  to 
have  lognog  slopes  of  about  -5/3.  Such  spectra  often  are  used  in  simulating  an  environment  or 
predicting  atmospheric  structure,  Fast  Fourier  transform  analysis  provides  a  means  for  filtering 
white  noise  with  spatial  fllLers  to  simulate  a  stationary  time  or  spatial  data  set.  In  many 
apphcations  the  transform  technique  provides  adequate  prr''''';ing  speed.  For  example, 
repeatedly  using  the  double  precision  IMSL  routine  DF2TCF,  ?  .  point  transform  averages 

6.91  ms  on  tlie  Phillips  Laboratory  Convex  model  210  computer,  v^cophysical  data  however,  are 
multidimensional,  and  oiteu  exceed  1024  points.  Three  dimensional  Fourier  simulation  of  atmos¬ 
pheric  structure  that  realistically  allows  for  variable  coherence  length  scales,  RMS  huctuauon 
levels,  and  spectral  slopes  would  require  months  of  execution  time  on  a  work  station. 

The  Phillips  Laboratcry  Strategic  High  Altitude  Atmospheric  Radiance  Code  (SHARC)’  uses 
first  principles  to  calculate  point -to-space  and  limb  \hewing  atmospheric  background  infrared 
radiance  and  transmittance.  Real  atmospheric  infrared  background  perturbations  occtu'  from 
fluctuations  in  temperature  and  density  of  the  contributing  molecular  species.  Version  4  of  the 
SHARC  code  emhsions  a  capability  to  evaluate  radiance  structure  from  estimated  variances  in  the 
standai'd  temperature  and  density  profiles.  To  provide  a  realistic  hut  practical  two-dimensional 
structure  scene  capability  will  require  creative,  efficient,  and  tested  algorithms.  The  purpose  of 
this  report  is  to  study  the  possibility  of  producing  s>mthetic  structure  from  autoregression 
analysis  as  contrasted  wth  the  common  Fourier  method  with  a  Mew  toward  reduemg  the 
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computational  burden.  Only  one-dimensional  analyiis  is  Seated  hire  but  a  subsequent  report  will 
present  a  two-dimensional  approach. 


a.  THEORY  OF  AUTOREGR^SIVE  DIGITAL  SPECTRAL  BSTralATION 

This  section  is  a  brief  review  of  autoregressive  digital  spectral  estimation  as  presented  by  S. 
Lawrence  Marple^  and  Steven  M.  Kay^. 

An  autoregressive  moving  average  (ARMA)  model  for  a  discrete  time  series,  x(n),  that 
approximates  deterministic  and  stochastic  processes  can  be  represented  by  die  filler  linear 
difference  equation: 

P  7 

x{n)  =  ■-ya{k)x{n-k)-i-'^b{k]z{n-k) 

i-=l  i=0 

in  which  xtn)  is  the  output  sequence  and  £{n)  a  white  noise  input  driving  sequence.  The  a(k)  and 
b(k]  coef&cients  form  the  autoregressive  and  moving  average  portions  of  the  ARMA  model 
respectively.  A  z  transform  analysis  of  the  difference  equation  shows  that  the  ARMA  power 
spectral  density  (PSD)  is: 

•Par.ua  (/)  “  Tpw 
where. 

Mf)  =  1  +  X  exp(-27:i/?cT) 

k=\ 

and, 

=  1  ^^b(fc)exp(-27:i/kr). 

Ic=l 

T  is  the  sampling  interval  and  is  the  variance  of  the  white  noise  process.  If  all  the  moving 
average  coefficients  are  zero,  except  b(0)  =  1.  then 

p 

x(n)  =  -^a(k:)x(n-  k)+  e(n) 
i<=i 

and  the  process  is  strictly  autoregressive  of  order  p.  The  autoregressive  PSD  becomes. 

^Marple.  S.L ,  1987  Digital  Spectral  Analysts  uiith  Applications.  Chapter  6.  Prentlce-Hal!.  Englewood  Cliffs.  New 
Jersey. 

^Kay.  Steven  M  .  1938  Modem  Spectral  Estlrration,  Theory  &  Application.  Prentice-Hall.  Englewood  Cliffs.  New 
Jersey. 
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Pmil)  -  «  T  Xr«(k)exp(-a^ykT) 


fes- 


whsre  rjg.  is  the  autocorrelation  sequence  at  k  and  Ow  is  the  variance  of  e(n).  The  relationship 
between  the  autocorrelation  sequence  and  the  pure  autoregressive  model  is 


Xa(ky^(m-k) 

ksi 

^a(k>-j^(--k)  +  p. 


k''l 


form>0 

for  m  =  0 
for  m<Q 


This  expression  may  be  evaluated  for  the  p+1  lag  indices  0  <  -•  <  p  by: 


■  A.(-P)'' 

"  1  ■ 

'P,' 

'•«(!)  ^xxiO)  ■' 

^^(1) 

0 

jM  r„{p-l)  • 

■■  At(O)  > 

,a(p)J 

This  expression  forms  the  autoregressive  Yule-Walker  equations.  Given  the  autocorrelation 
sequence  for  lags  0  to  p,  the  autoregressive  coefficients  may  be  found  from  the  above.  Since 
-k)  =  r'jQ-fk),  the  autocorrelation  matrix  is  both  Toeplitz  and  Hermitian.  A  standard  "Levinson" 
algorithm,  which  takes  act/antage  of  the  Hermitian-Toeplitz  matrix  equation,  was  employed  to 
solve  for  the  AR  parameters.  The  same  Yule-Walker  equations  also  occur  if  we  attempt  to  solve  the 
problem:  find  the  ''best",  in  a  least  square  sense,  set  of  equations  that  determine  the  coefficients 

p  - - - - 

Uiat  predict  Sl  from  x{n)  =  -^a{k)xin-k).  where  (.x(n)-  .tfn))"  =  p.,.  In  summary,  the  Levinson 

v=i 

recursion  computes  sets  of  coefficients  {ai(l).  pj).  (02(1).  a2(2),  P2! .  (apd),  ap(2) .  ap(p),  pp) 

where  the  final  set  at  order  p  is  Uie  desired  solution  of  the  Yule-Walker  expressions.  For  the  AR(p) 
process.  ap(T  =a(i)  for  i  =  1,  2,  3 . p  and  pp  is  the  minimum  prediction  error,  that  is,  p.j^r  =  pp  = 

pjTiin  ^fx''[n](x[n]- x[n.])j.  The  algorithm  is  initialized  by; 


Ax[0] 


«:[!]  = 


with  the  recursion  for  k  =  2,  3 . p  gwen  by 
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Pi-1 


Ok  [i]  =  ajfe-i  [i]  +  %[^]flM  [k-i]  !  -  1. 2 . ^'  - 1 

Pi  “  (^~!'^i[^'ll  )Pi-l 

The  ajjikl  coefficients  are  known  as  reflection  coeffidents. 


3.  AFPLICATION  TO  ATMOSPHERIC  POWER  SPECTRAL  DENSITIES 

Atmospheric  power  spectral  density,*  functions  often  are  modeled  by  power  law  functions  of  the 
form‘''^ 


PSD(k)^ 


The  constant  C  is  found  from  the  definition  that  the  total  variance  of  the  associated  time  series, 
ct^,  is  equal  to  the  area  under  the  PSD®. 


a^  =  2j?5D(l:)dJi  =  2cj‘ 


r  dk 


Vnr(v) 

\  2 


•:s2C- 


so  that  the  PSD  model  is 


PSD{k)  = 


I  2vp^  1  ''i 

a  a  Pi  v-i-- 
V  2  j 


i/rtr{\-){a~  +k‘') 


ZV*'/! 


The  relationship  betiveen  the  frequency  domain  PSD  and  the  time  or  spatial  domain 
autocorrelation  function  is  specified  by  their  Fourier  transform  pans.  Thus  the  autocorrelation 
function.  (ACF).  for  the  real  even  PSD  function  is 


^atarskl.  V.I..  196i  WcfC  Propagation  in  a  Turbulent  Afedtiun.  McGraw-Hill.  , 

^Putterman,  W.i  .  SchwelUer.  E.L..  and  Newt.  J.E.  1991  "Estimation  of  Scene  Correlation  Lengths". 
Chaxact''rlzadoij,  Propagation,  and  SLirulatlon  of  Sources  and  Backgrotinds.  Proceedings  SPIE  ■  The  International  Society 
of  Optical  Englneei'lng,  V  1486.  pp  127- 140.  Orlando.  Florida. 

®For  the  integral,  see  for  example.  GradshlejTi.  l.S.  and  Kyzhlk,  I.M..  1965  Table  of  Inlegrals  Senes  and  Products,  eq 
3.241.4.  Academic  Press 
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o^rf  V +|- V'’  Qot{2mk)dk 
ACFix)mIFmPSD)^i[  ■ — -4-  -  ‘i/- 


The  intepaJ  is  evaluated  using  the  Bessel  function  of  the  second  kind  of  fractional  order.  to 
be' 


ACFC-O  i 


r(v)' 


Ejqpressions  for  two  and  Uii'ee  dimensional  power  spectral  densities  can  be  evaluated  from  the 
model  autocorrelation  function  by  performing  two  and  dnee  dimensional  Fourier  transforms  of 
ACF(x).  For  example,  over  the  valid  range  for  v,  the  two-dimensional  isotropic  power  spectral  den¬ 
sity  is 


F(t-)- jy^j|cosft“  f)ACF(?'icP^ 


v+i 


r.{a^+ftT 

while  the  three-dimensional  isotropic  power  spectral  density  is 
<!>(?)  =  -^^^|JJcos(£'r)ACF(r>fr.  .  . 

n^T{v){a^  ^ 

The  parameter  "a"  can  be  expressed  in  terms  of  the  mtegral  scale  or  equivalent  width  of  the  ACF. 
An  equivalent  width  is  defined  as  the  area  of  a  function  dmded  by  its  central  ordinate®,  or 

'  7(0)““ 

Thus  the  large  scale  correlation  length  is  defined  by  integrating  ACF(x)  over  all  positive  values 
of  X. 


Ir  = 


f  ACF(x)alx 

Jo 


a" 


■^Ibld..  Eq  £  -'32.5 

®Bracewell.  RN.,  1978  Ttie  Fourier  Transjortn  and  Its  AppUcaitans.  McGraw-HJl.  New  York 


-wntre  autocorrelation  at  zero  separatioo.  Obieiving  that  the  equivalent  width  of  a 

function  is  equal  to  the  inverse  of  the  equivalent  width  of  its  transform®,  the  correlation  lengUi 
may  le  written  as 

L  P^g(O)  ■■■  PSDjO) 

2j^PSDik)dk 

so  that, 


4 


°  V  ir 

2cT^-s/Rr(v)a^''-*' 


r 

n  V 


2Vnr(v)fl 


or 


a  = 


2V?r(v 

rfv*i' 


These  relationships  are  visually  compared  in  Figure  1.  The  upper  left  quadrant  of  the  graph  shows 
a  log-log  plot  of  the  PSD  function  for  a  slope  of  -3.  ~  0.2.  =  5.6  and  a  =  0.05.  The  upper  right 
quadrant  shows  a  linear  plot  of  the  PSD  with  L(,  and  calculated  from  the  curve.  The 
autocorrelation  function  is  plotted  in  the  lower  left  quadrant  and  also  shows  and  calculated 
from  the  ciuve.  Quite  obviously  the  parameter  "a"  depends  upon  the  %'alue  of  v,  such  that  "a" 
coincides  with  the  roll  over"  or  "comer  frequency"  of  the  PSD.  The  correlation  length,  on  the  other 
hand,  corresponds  to  a  much  higher  frequency  (and  smaller  PSD).  As  shoi\Ti  in  the  lower  right 
quadrant  of  Figure  1.  the  parameter  "a"  assumes  a  value  near  the  frequency  of  the  peak  power 
while  l/Lj  assumes  a  value  on  the  far  tail  of  the  power  curve, 


4.  IMPLEMENTATION 


In  practice  a  given  power  law  PSD  was  constructed  from  the  PSD  and  ACF  models; 

2alAx 


PSD{f]  = 


c  a  r.v.-J 


VS'r(v)(o'-^/’) 


i  1 


ACF{x)  = 


K,.{2ru3x] 

r(v) 


^Ibld 
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where, 

cj  ■  residual  error 

^  ■  spacing 

aj[  ■  model  parameters 

f  H  frequeiiQf  (Hz) 


Then,  a  spatial  data  sequence  was  simulated  using  a  normally  dlstrtbuted  set  of 
pseudorandom  numbers  G(J).  The  simulation  of  M  sequence  values  proceeded  as  follows; 

SetY(l)  =  Y(2)  =  Y(3)  ^  ■  =Y(N)  «  0.  Then  for  J  =  N+i.  N’+2.  M+IOOO.  create  Y(J)  from 

Y{J)  ^  G{J)a^-'y  a^Y{J-i).  The  simulated  sequerne  is  then  contained  in  Y(lOOl).  Y(1002),  ■ 

(=1 

Y(M+1000), 

After  simulating  such  a  discrete  spatial  sequence  of  a  given  PSD,  the  PSD  of  the  simulated 
series,  was  .checked  as  follows.  Using  a  forward  and  backward  estimation  method^®,  a  set  of  NN  bi 
coefficients  were  found  to  estimate  the  simulated  PSD  from  the  formula. 


P5D(/)- 


-26^  Ax 


SN 


(s! 


where  Ax  is  the  same  spacing  used  to  generate  the  synthetic  series. 

The  forward  and  backward  metliod  soWes  the  following  least  squares  problem.  Given  M 
discrete  series  values  Y(.J).  for  J  =  1.2.  ■  .  M.  we  fmd  bj  values  that  minimize  ERR,  where, 

f  ^  VI  S  V''. 

ERR  =  \  '  Y{J)-^btY[J-i)\  j+i  ^  !  ^(7)- + i) 

■^/=.n7,'+'A  i=;  y  V  ^  J 


CTp  is  then  equal  to 


ERR 


2(iVf-A7V) 


5.  RESULTS 


^^Haykln.  Simon,  ed,,  1991  Advances  (;i  Spectrum  Analysts  and  Array  Processing.  Vol.  1.  Prentice  Hall.  Englewood 
CUfl's  New  Jersey.  .  pp.  155- 156. 
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The  following  discuision  Is  atoad  at  producing  a  practical  autoregresiivi  model  of 
atmospheric  structure  consistent  with  having  a  pra-asil^ad  power  spectral  densl^  and 
autocorrelation  function.  Correlation  lengths  and  d®  variances  are  taken  from  Strugala.  et.al." 

Figure  2  shows  an  8000  point  sequence  of  Gaussian  "white  noise"  having  a  spacing  of  100  m. 
mean  of  zero,  and  standard  deviation  of  O.OS5.  Figure  3  presents  a  histogram  of  tiie  Gaussian 
data  and  Figure  4  shows  the  "flat"  white  noise  spectra  and  associated  autocorrelation  function. 
The  theoretical  and  calculated  statistics  agree,  with  the  calculated  having  an  average  of  0.000117 
and  0  of  0.0B49,  Correlated  sequences,  histograms,  PSD's  and  ACF's  of  subsequent  plots  may  be 
compared  to  this  "unfiltered  sequence." 


Figure  5  is  typical  of  all  the  plots  in  this  report  showing  log-log  PSD's  in  tlie  left  panel  and 
ACF's  in  the  right  panel.  This  and  subsequent  plots  have  PSD's  measured  in 


(o  Temperature  /  Temperature)^ 
Wavenumber 


and  wavenumber  measured  in  km'b  Three  curves  usually  appear  in 


each  panel.  The  solid  unmarked  curves  are  the  "theoretical"  or  inputted  PSD  or  ACF.  The  curv'es 
marked  by  an  X  are  tire  autoregressive  predicted  PSD's  and  ACF's.  and  the  curves  marked  by  an 
open  square  are  PSD's  and  ACF's  derived  from  the  simulated  data  sequences.  Except  where 
otherwise  noted,  the  data  spacing  is  100  m  and  the  Nyquist  frequency  is  5  km*’.  The  input 
parameters  for  Figure  5  are  1.'75  km.  »  0.00102.  spectral  slope  (S)  ^  -B/S.  Using  six  linear 
predictor  coefficients,  the  predicted  values  are  1.44  km  from  tlie  PSD  and  1.50  km  from  the 
ACF.  and  d2  =  0.00107  from  the  area  of  the  PSD.  Using  12  coefficients  to  determine  the  simulated 
results,  the  values  for  the  simulated  curves  are  Lj.*  1.38  km  from  the  PSD  and  1.44  km  from  the 
ACF,  and  =  0.00106  from  the  area  of  the  PSD.  Visually  examming  Figure  5,  clearly  shows  that 
usmg  6  coefficients  for  the  predictor  and  12  coefficients  for  the  simulated  values  reproduces  the 
theoretical  slope  quite  well  except  at  the  highest  wavenumbers.  Also,  at  low  wavenumbers,  the 
predicted  and  simulated  PSD's  show  a  “16  percent  DC  bias  from  the  theoretical.  Part  of  the  hlgli 
wavenumber  divergence  Is  due  to  the  sharpness  of  the  autocorrelation  function  at  zero  lag.  In 
fact,  for  the  specified  input  spectral  slope  of  “5/3.  a  sharp  cusp  exists  at  zero  lag.  This  cusp  intro¬ 
duces  an  error  into  the  estimate  of  the  linear  prediction  coefficients  tliat  affect  the  PSD  at  high 
wavenumbers.  To  lest  this  h\q)Othesis,  the  first  value  of  the  autocorrelation  function  at  zero  lag 
was  modified  by  setting  it  equal  to  a  \'alue  linearly  extrapolated  from  the  second  and  third  values. 
This  modification  led  to  the  curves  in  Figure  6.  Here  we  see  much  better  agreement  between  the 
theoretical  and  modeled  curves  at  high  wavenumbers  and  improvement  at  low  wavenumbers.  Im¬ 
provement  in  the  spectral  parameters  is  also  e^ndent.  Lj.  improves  to  1.54  km  for  the  predicted 
PSD  and  1.60  km  for  the  predicted  ACF.  Figure  7  shows  a  sequence  of  8000  data  points  simulated 
over  800  km.  When  compared  witlr  Figure  2.  the  -5/3  spectral  filtering  is  obvious.  To  see  if  the 
simulated  data  retains  the  Gaussian  PDF,  a  histogram  of  the  simulated  data  is  examined  in 
Figure  8.  The  figure  shows  that  the  variance  of  the  simulated  data  =  0.00102.  matching  the  input 
parameter,  o^.  Also,  the  histogram  closely  follows  the  shape  of  the  theoretical  PDF  as  shown  by 
the  solid  cun-e.  The  effects  of  using  more  or  fewer  predictor  coefficients  are  shown  in  Figures  9 


^  ^strugala.  Lj^..  Nev.!,  J.E..  Fuctenr.an.  W.I .  Schweitzer.  E.L..  Herman.  B.J..  and  Sears,  RD.  1991  Development  of 
High  Resolution  Statistically  Non-statlonaiy  Infrared  Earthllmb  Radiance  Scenes.  Charactertzation.  Propagation,  and 
Sirnulatlon  of  Sources  and  Backgrounds,  Proceedings  SPIE  -  The  International  Soclew  for  Optical  Engineering.  VI 486.  np 
176-187.  Orlando.  Florida. 
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and  10.  figure  9  is  calculated  for  only  one  predictor  coefficient  while  figure  10  is  calculated  using 
20  predictor  coefficients.  Except  for  the  bias  at  low  wavenumbert,  using  only  one  coefficient  for 
the  predictor  does  remarKably  well  for  the  specffied  input  parameters.  However,  with  tlie  use  of  20 
coefficients,  any  dlffeience  between  the  theoretical,  predicted,  and  slmirlated  curves  is  practically 
erased.  Resultant  values  for  die  predicted  and  simulated  spectral  parameters  also  improve.  For 
example,  Lg  improves  to  1.63  km  for  the  simulated  PSD  and  1.70  km  for  the  simulated  ACF. 
Figure  11  provides  a  visual  perspective  of  how  the  20  corresponding  reflection  coefficients 
decrease  as  a  function  of  reflection  coefficient  number,  Ordinarily,  the  reflection  coefficient 
decreases  smoothly  and  monotonically,  However,  introducing  the  modification  to  tlie  ACF  that 
was  mentioned  above  causes  the  second  reflection  coefficient  to  dip  below  the  smooth  value, 
Evidently  this  helps  compensate  for  the  cusp  in  the  ACF, 

One  may  see  the  effects  of  increasing  the  value  of  the  correlation  length  Lg,  in  Figures  12-14. 
These  cuives  are  calculated  for  Lc  =  10  km.  =  0.0049,  and  S  =  -5/3.  Using  6  coefficients  for 
the  predicted  PSD  and  12  coefficients  for  the  simulated  PSD.  produces  the  curves  of  Figures  12 
and  13.  Figure  12  illustrates  the  increased  difficulty  in  simulating  a  PSD  and  ACF  ha\ing  a  large 
Lg  with  a  small  number  of  coefficients.  Figure  13  clearly  shows  the  effect  of  the  larger  correlation 
length  on  the  simulated  spatial  sequence.  As  expected,  greater  srnoothmg  (data  coirelatlon) 
results  from  using  the  larger  scale  size.  Figure  14  Illustrates  the  improv'ement  in  the  predicted 
and  simulated  PSD  and  ACF  by  using  20  predictor  coefficients  for  both  the  predicted  and 
simulated  curves. 

F'igures  15-20  show  the  result  of  repeatmg  the  calculations  for  a  spectral  slope  of  -2.  Figures 
15-17  have  inputted  spectral  parameters  of  Lg  =  1.75  and  =  0.00102  while  Figures  18-20  have 
Lg  a  10  and  cj2  «  0.0049.  For  a  slope  of  -2,  using  6  coefficients  for  the  predicted  PSD  and  12 
coefficients  for  the  simulated  PSD.  gives  good  agreement  between  theoretical,  predicted,  and 
simulated  curves,  even  for  Lc  =  10.  For  a  theoretical  Lc=  1.75.  ^  0,00102.  the  predicted  values 

are  1.73  from  the  PSD  and  1.75  from  the  ACF.  and  =r  0.00104  from  the  area  of  the  PSD, 
The  values  for  the  simulated  curves  are  Lc=  1.65  from  the  PSD  and  1,67  from  the  ACF,  and  = 
0.00103  frjm  the  area  of  the  PSD.  Figure  16  shows  the  sequence  of  8000  data  points  simulated 
over  800  km  for  the  -2  spectral  slope  and  Lc=  1.75.  When  this  result  is  compared  with  Figure  7, 
the  increased  filtering  is  evident.  Again  the  simulated  data  retains  the  Gaussian  PDF  as  shown  in 
the  histogram  of  Figure  17.  The  figure  shove’s  that  the  variance  of  the  simulated  data  =  0.00102. 
matching  the  input  parameter,  o^.  As  shown  in  Figure  18.  inputting  a  theoretical  Lc=  10  and  = 
0.0049  gives  the  predicted  values  of  Lc=  9.98  from  the  PSD  and  10  from  the  ACF,  and  = 

0.00491  from  the  area  of  the  PSD.  The  values  for  the  simulated  curv'es  are  Lc=  9.37  from  the  PSD 
and  9.39  from  the  ACF.  and  =  0.00472  from  the  area  of  the  PSD.  Figure  19  shows  the  se¬ 
quence  of  8000  data  potn.ts  simulated  over  800  km  for  tlie  -2  spectral  slope  and  Lc=  10.  Figure  20 
again  show's  that  the  histogram  of  the  simulated  data  retains  the  Gaussian  PDF  and  that  the 
v'ariance  of  the  simulated  data  =  0.00468,  approximating  the  input  parameter, 

Figure  21  repeats  the  calculations  for  a  theoretical  sperh-aJ  slope  of  -3.  Lc=  1.75,  and  = 
0.00102.  Up  to  this  point  the  autocorrelation  function  has  been  modified  at  zero  lag  as  discussed 
above.  However,  as  evident  in  Figure  21.  for  slopes  of  -2  and  greater,  the  modified  ACF  results  in 
a  divergence  of  the  predicted  PSD  at  high  wavenumbers  whereas  the  unmodified  ACF  produces 
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good  fldellt*/  at  high  wavenumbers.  For  this  reason,  use  of  the  modlfted  ACF  is  appropriate  only 
for  tile  -6/3  spectral  slope.  For  slopes  of  -2  and  larger,  subsequent  analysis  will  use  the 
'onmodified  ACF  exclusively. 

Figures  22-27  show  the  result  of  repeating  the  calculations  for  a  spectral  slope  of  -3.  Figures 
22-24  have  inputted  spectral  parameters  of  Lg  ^  1.75  and  o'^  =  0.00102  while  Figures  26-27  have 
Lg  10  and  0^  =  0.0049.  For  a  slope  of  -3.  using  6  coefficients  for  the  predicted  PSD  and  12 
coefficients  for  the  simulated  PSD.  good  agreement  obtains  between  theoretical,  predicted,  and 
simulated  curves,  although  tliere  is  less  agreement  at  low  wavenumbers  for  Lg  *  10.  For  a 
tlieoretical  L^s  1.75,  =  0.00102.  the  predicted  values  are  Lc=  1.88  from  the  PSD  and  1.88  from 

the  ACF,  and  ^  0.00102  from  the  area  of  the  PSD.  The  values  for  the  simulated  curves  are  Lgs 
1.79  from  the  PSD  and  1.79  from  the  ACF,  and  -  0.00102  from  the  area  of  the  PSD.  Figure  23 
shows  the  sequence  of  8000  data  points  simulated  over  800  km  for  the  -3  spectral  slope  and  Lg= 
1.75.  Comparisons  with  Figure  16  show  the  -3  versus  -2  increased  filtering.  Again  the  simulated 
data  retains  the  Gaussian  PDF  as  shown  in  the  histogram  of  Figure  24.  The  figure  shows  that  the 
variance  of  the  simulated  data  =  0,00102,  matching  the  input  parameter,  o-^.  As  shown  in  Figure 
25,  inputting  a  theoreUcal  L(,=  10  and  =  0.0049,  gives  predicted  values  of  6^=  16.6  from  the 
PSD  and  16.6  from  the  ACF,  and  =  0.00490  from  the  area  of  the  PSD.  The  values  for  the 
simulated  curves  are  6^=  16.1  from  the  PSD  and  16.1  from  the  ACF,  and  -  0.00473  from  the 
area  of  the  PSD.  Figure  26  shows  the  sequence  of  8000  data  points  simulated  over  800  km  for  the 
-3  spectral  slope  and  10.  Figure  27  again  shows  that  the  histogram  of  the  simulated  data 
retains  the  Gaussian  PDF  and  that  the  variance  of  the  simulated  data  =  0.00471,  approximating 
the  input  parameter,  r'igure  28  shows  the  improvement  that  is  made  by  using  12  predictor 
coefficients.  In  this  case,  the  predicted  values  are  13  from  the  PSD  and  13  from  the  ACF.  and 
o2  =  0.00490. 

The  final  spectral  slope  examined  in  this  report  is  for  a  slope  =  -4.  Figures  29-34  present  the 
results  of  repeating  the  calculations  for  a  S  =  -4.  Figures  29-31  have  inputted  spectral  parameters 
of  Lc  -  1.75  and  =  0.00102  while  Figures  32-34  have  =  10  and  =  0.0049.  In  tliis  case, 
using  6  coefficients  for  the  predicted  PSD  and  12  coefficients  for  tlie  simulated  PSD  ynelds  good 
agreement  between  theoretical,  predicted,  and  simulated  curves.  For  a  theoretical  1.75,  = 

0.00102.  the  predicted  values  are  1.75  from  the  PSD  and  1.75  from  the  ACF,  and  = 
0.00102  from  the  area  of  the  PSD.  The  values  for  the  simulated  curves  are  L(;=  1.67  from  the  PSD 
and  1.67  from  the  ACF,  and  ~  0.00102  from  the  area  of  the  PSD.  Figure  30  shows  the 
sequence  of  8000  data  points  simulated  over  800  km  for  the  -4  spectral  slope  and  1.75. 
Comparisons  with  Figure  23  show  the  -4  versus  -3  increased  filtering  effect.  Again  the  simulated 
data  retains  the  Gaussian  PDF  as  showm  in  the  histogram  of  Figure  31.  The  figure  shows  that  the 
variance  of  the  simulated  data  =  0.00102,  matching  the  input  parameter,  as  shown  in  Figure 
32,  inputting  a  theoretical  10  and  =  0.0049,  gives  predicted  values  of  Lq-  9.98  from  the 
PSD  and  9.98  from  the  ACF.  and  =  0.00490  from  the  area  of  the  PSD.  The  values  for  the 
simulated  curves  are  Lc=  10.1  from  the  PSD  and  10.1  from  the  ACF.  and  =  0.00463  from  the 
area  of  the  PSD.  Figure  33  shows  the  sequence  of  8000  data  points  simulated  over  800  km  for  the 
-4  spectral  slope  and  10.  Figure  34  again  shows  that  the  histogram  of  the  simulated  data 
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retains  the  Gaussian  PDF  and  that  the  variance  of  the  simulated  data  *  0.00461,  appro^dmating 
the  input  parameter, 

Up  to  this  point  the  anatysis  has  been  based  on  a  data  spacing  of  100  m,  The  effect  of 
doubling  Oie  spacing  Is  demonstrated  in  Figures  35  and  36,  These  plots  have  inputted  spectral 
parameters  of  S  -  -5/3  and  ^  B.  Figure  35,  which  was  calculated  for  a  data  spacing  of  200  m. 
may  be  compared  with  Figure  36  which  has  a  data  spacing  of  100  m.  Of  course,  doubling  the  data 
spacing  from  100  to  200  m  halves  the  spatial  resolution  and  reduces  the  Nyquist  frequency  from 
5  tai'*  to  2.b  knr^.  If  reducing  resolution  is  acceptable,  the  positive  effect  is  the  improvement  in 
the  predicted  correlation  length,  For  example,  by  using  6  linear  predictor  coefficients,  the 
predicted  Lc=  3.6  km  for  a  spacing  of  100  m  but  4.28  km  for  a  spacing  of  200  m. 

6.  CONCLUSION 

Geophysical  phenomena  within  a  specified  domain  often  are  chmacterized  by  smooth 
continuous  power  spectral  densities  having  a  negative  power  law  slope.  The  association  of  one-, 
two-,  and  three-dimenslonaJ  geophysical  spectral  densities  with  a  given  autocorrelation  function 
was  reviewed  and  the  autoregressive  methods  of  modem  spectral  estimation  were  explored.  An 
autoregressive  alternative  to  the  more  common  Fourier  transform  analysis  was  studied  with  a 
goal  of  reducing  the  enormous  computational  burden  of  generating  synthetic  structure  scenes 
from  Gaussian  random  number  sets.  A  6  coefficient  AR  run  to  simulate  1300  points,  using  the 
Phillips  Laboratory  model  210  Convex  computer,  showed  an  average  execution  time  of  0.462  ms. 
an  imprm'ement  of  a  factor  of  6.91/0.462  over  the  fast  Fourier  transform.  In  generating  two- 
dimensional  synthetic  arrays,  greater  savings  would  result  due  to  the  8-fold  symmetry  of  the 
quarter-plane  AR  coefficients.  It  was  demonstrated  how  the  resolution  and  accuracy  of  predicted 
and  simulated  data,  their  PSD's  and  ACF's,  and  their  parameters,  change  with  spectral  slope, 
correlation  length,  data  spacing,  and  prediction  order.  In  particular,  six  linear  prediction  coeffi¬ 
cients  are  all  that  is  necessary  in  many  cases  to  generate  a  synthetic  spatial  sequence  that 
retains  a  specified  power  law,  correlation  scale,  variance,  and  probabiUt>’  distribution  function. 
For  a  desired  spectral  slope  of  -5/3  one  should  employ  an  autocorrelation  function  having  a 
modified  zero  lag  value.  For  an  equeil  linear  predictor  order,  structured  data  havmg  smaller 
correlation  length,  larger  spectral  slope,  or  reduced  data  resolution,  have  greater  fidelity  to 
specified  characteristics  than  those  generated  for  small  slope,  larger  correlation  length,  or  higher 
resolution.  However,  fidelity  often  can  be  reestablished  by  increasing  the  spectral  order  from  6 
coefficients  to  12  or  in  the  worst  case  20  coefficients,  Sir.ce  the  maximum  vertical  correlation 
length  over  relevant  SHARC  altitudes  is  reported  to  be  approximately  10  km,  a  vertical  resolution 
of  100  m  can  be  achieved  with  a  minimum  of  coefficients.  However,  smce  maximum  horizontal 
correlation  lengths  are  reported  to  be  approximately  85  km,  horizontal  resolution  must  be 
sacrificed  or  else  many  more  coefficients  must  be  used  to  achieve  fidelity  at  both  low  and  high 
frequencies.  Where  low  frequency  simulation  is  unimportant,  one  may  retreat  to  a  minimum 
number  of  coefficients  and  still  achieve  fidelity'  at  mid  and  high  frequencies. 


Autoregressive  (aR)  and  autoregressive  moving  average  (ARMA)  modeling  should  be 
considered  in  creating  large  atmospheric  structured  scenes.  A  subsequent  report  will  address  the 
potential  for  using  two-dimensional  ARMA  modeling  to  generate  two  and  tlirce-dimensional  struc¬ 
tured  scenes. 
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.  Sample  Gaussian  riindom  number  sequence,  mean  =  0..  standard  deviatloo  =  0.055,  spacing  100  m 
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Figure  9.  Left  panel,  theoretical  PSD  input  (unmarkedK  autore^esstve  predictor  PSD  (marked  by  xj;.  Right  panel,  corresponding 
autocorrelatfon  functions.  =  1.75  km.  S  =  5/3  a  =  1.O2E-03,  spacing  =  100  m.  Modified  theoretical 
autocorrelation  function  at  lag  =  0.  One  predictor  coefficient 
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spacing  =  100  m.  Theoretical  autocorrelation  function  modified  at  lag  =  0.  Tiventy  predictor  coefflcients 
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F'lfinre  13.  Autoregressive  simulated  structure  sequence.  —  10  km,  S  —  -5/3.  a  -  4.'BE-03 
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Figure  14.  Left  panel,  theoretical  PSD  Input  (unmarked^  autoregressive  predictor  PSD  (marked  by  x}.  and  ^mulated  PSD  {marked 
by  small  square).  Right  panel,  corresponding  autocorrelation  functions.  L^,  -  10  km.  S  =  -5/3.  a  =  4.9E-03.  spacing  = 
100  m.  Theoretical  autocorrelation  function  modified  at  lag  =  0.  Tlwenty  predictor  coefficients 
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Figure  15.  Left  panel,  theoretical  PSD  Input  (unmarked},  autoregressive  predictor  PSD  (marked  by  x},  and  |lniulated  PSD  (marked 
by  small  square].  Right  panel,  corresponding  autocorrelation  functions.  Lj.  =  1.75  km,  S  =  -2,  a  =  1-02E-03,  spacing 
100  m.  Theoretical  autocorrelation  function  modified  at  lag  =  O.  Six  predictor  coefficients 
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Figure  16.  Autoregressive  simulated  structure  sequence.  =  1.75  km,  S  =  -2.  a  =  I.02E  03 


Figure  17.  Histogram  of  autoregressive  simulated  structure  sequence.  L(.  =  1.75  km,  S  =  -2.  a  =  1.02E-03 
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Figure  19.  Autoregressive  siimilated  structure  sequence.  =  lO  km.  S  =  -2,  a  =  4.9E-03 


Figure  20.  Histogram  of  autoregressive  simulated  stnicture  sequence.  =  10  km,  S  =  2,  a  -  4.9E-03 
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Figure  22.  Left  panel,  theoretical  PSD  Input  (unmarked),  autoregressive  predictor  PSD  (marked  by  x),  and  |{mul'ated  PSD  (raairked 
by  small  square).  Right  panel,  corresponding  autocorrelation  functions.  L^.  =  1.75  km,  S  =  -3.  u  =  I.02E:-0-3.  spacing  = 
100  m.  Unmodified  theoretical  autocorrelation  function.  Six  predictor  coefficients 


Fli'ure  23.  Autoregressive  simulated  structure  sequence.  L^.  =  1.75-  km.  S  =  -3.  a  =  I  02E-03 
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Rgure  24.  Histogram  of  autoregressive  simulated  structure  sequence.  L^.,  =  1.75  km,  S  =  -3,  0  =  i.02E-03 
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Figure  25.  Left  panel,  theoretical  PSD  Input  (unmarkedK  autoregressive  predictor  PSD  (marked  by  xj.  an^  simulated  PSD  {marked 
by  small  square).  Right  pemel,  corresponding  autocorrelaLton  functions.  L^.  -  10  km.  S  =  -3,  o  =  4.9E-03.  spacing  = 
100  m.  Unmodified  theoretical  autocorrelation  function.  Six  predictor  coefficients 


UJ  liJ 


Q 

O 

O 

O 


\A  . 

2 

P 


£ 

—  o 
<J  X 


Si/. 

‘•a) 


X  a 
8  £ 
B  ui  ; 


G  C  S 


V.  2  C  X 


38 


Fljfure  26.  Autoregressive  simulated  structure  sequence.  Lj.  =  10  km,  S  =  -3;  o  =  4.9E-03 
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Fiflure  27.  Histogram  of  autoregressive  simulated  structure  sequence.  =  10  km,  S  =  -3,  tr  =  4.9E-03 
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Figure  2^.  Left  panel,  theoretical  PSD  input  (unmarked),  autoregressive  predictor  PSD  (marked  by  x),  and  |tmulated  PSD  Cnnarked 
by  small  square).  Right  panel,  corresponding  autocorrelation  functions.  =  1.75  kin,  S  =  -4,  a  =  I.02E:-03„  spacing  = 
100  m.  Unmodified  theoretical  autocorrelation  function.  Six  predictor  coefficients 


Figure  30.  Autoregressive  simulated  structure  sequence.  -  1.75  km.  S  =  -4.  c  =  1.02E-03 


Figure  3 1 .  Histogram  of  autoregressive  simulated  structure  sequence.  Lc  =  1.75  km,  S  =  -4.  o  =  L02E-03 
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Autoregressive  simulated  structure  sequence.  =  10  km.  S  =  -4.  o  =  4.9E-03 
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Figure  35.  Left  panel,  theoretlral  PSD  Input  (unmarked^,  autoregressive  predictor  PSD  (marked  by  xK  Right  panel,  corresponding 
autocorrelation  functions.  Lp  =  5  km,  S  =  -5/3,  spacing  =  200  m.  Tlieorellcal  autocorrelation  function  modified  at  Iag  = 
0.  Six  predictor  coefTlctents 
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Figure  36.  Left  panel,  theoretical  PSD  input  (unmarkect).  autoregressive  predictor  PSD  (marked  by  x}.  Right  panel  correspondfog: 

autocorrelation  functions.  =  5  km,  S  =  -5/3.  spacing  =  HX>  m.  Tlieoretical  autocorrelation  function  modifled  at  I’ag  = 
0.  Sl.r  predictor  coefTlcients 
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Appendix 

Extensions  of  AR  Methods  to  Non>Stationa»  y 
Data  and  to  Two  and  Three  Dimensions 


Ai.  nonjStationary  data 


The  methods  employed  here  easily  can  be  extended  to  simulate  a  non*stationaiy  discrete 
spatial  series.  If  Lg,  or  spectral  slope  are  c.Uowed  to  vary  slowly,  sets  of  af  coefficients  and 
associated  values  for  each  point  k  can  be  evaluated,  where  Is  the  1^  coefficient  ar  point  k 
with  the  simulated  results  starting  at  k  *=  1.  Y(k)  can  then  be  simulated  by  setting  Y(-L)  «  Y(*I>1)  - 
Y(L+2)  a  ..  a  Y(-L+N-l)  s  0  (where  L  Is  some  number,  say  100,  and  N  is  the  number  of 
coefficients).  With  k  *  -L+N,  -L+N+l,  •  l,  the  iterative  expression  for  Y(k)  is: 

.V 

Y(k)-u,G(fc)-5^a<‘Y(k-.U 

(si 


and  for  k  >  1,  Y(k)  is: 


a2.  two  and  three  dimensional  simulations 

The  methods  employed  in  the  text  can  be  e>d:ended  to  simulate  horizontal  and  vertlcail  data, 
Assume  we  wish  to  simulate  data  ai  altitudes  k  s  i,  2.  ■■■.  Then  using  either  a  Fourier  transform 
method  or  the  linear  prediction  technique,  one  may  form  a  simulated  sequence  whose  horizontal 
PSDs  are  the  PSD  from  the  L^.  and  slope  at  each  horizontal  altitude  with  a  set  to  1.  In  addition, 
one  must  form  simulated  values  at  <  Itltudes  -L,  -Dfl.  •••,  0  from  and  slope  (horizontal)  for  k  *1 
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and  for  o  set  to  1.  Let  us  call  the  reJultlng  v^ues  2(k-j)  where  k  is  the  altitude  and  J  is  a  point 
number  0  could  represent  a  point  on  a  horizontal  line  or,  by  eirtenslon,  a  point  in  a  two» 
dimensional  horizontal  sheet).  TTien  the  final  data  having  the  desired  horizontal  and  vertical 
PSD's.  F(kj),  is  formed  by  filtering  Z{kJ)  as  follows. 

Porks  -L.  ^L+1,  -L+2.  1 

S’ 

F(k.j)  =  o,Z(lc.j)-£a]Z(k-i,j) 

(si 

where  CT]  and  a}  are  the  values  found  using  tlie  Levinson  algorithm  with  rjgj  detennined  from  a 
and  the  vertical  values  of  L^  and  slope  at  altitude  ==1. 

Then  for  k  >  1, 


F(k.  j)  =  a^Zik.J)  -  £a|^Z(k  -  i.  j) 

(=1 

where  and  are  the  values  found  using  tire  Levinson  algorithm  with  r^x  determined  from  a 
and  the  vertical  values  of  L^.  and  slope  at  altitude  k. 

In  the  above  analysis,  it  should  be  noted  that  with  careful  programming  only  N+i  altitudes 
need  to  be  stored  at  any  time  during  the  computations.  This  can  reduce  the  necessary  storage  by 
a  factor  of  10  or  more. 
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